# cloglogsimcf QoI calculation code
# adapted from simcf package by Chris Adolph
# Inhwan Ko, Mar 8 2022

cloglogsimev <- function(x,b,ci=0.95,constant=1) {
  if (any(class(x)=="counterfactual")&&!is.null(x$model)) {
    x <- model.matrix(x$model,x$x)
  } else {
    if (any(class(x)=="list")) x <- x$x
    if (is.data.frame(x)) x <- as.matrix(x)
    if (!is.matrix(x)) {
      if (is.matrix(b)) {
        x <- t(x)
        if (!is.na(constant)) {
          x <- append(x,1,constant-1)
        }
      } else {
        x <- as.matrix(x)
        if (!is.na(constant)) {
          x <- appendmatrix(x,rep(1,nrow(x)),constant)
        }
      }
    } else {
      if (!is.na(constant)) {
        x <- appendmatrix(x,rep(1,nrow(x)),constant)
      }
    }
  }
  
  esims <- nrow(as.matrix(b))
  
  nscen <- nrow(x)
  nci <- length(ci)
  res <- list(pe = rep(NA, nscen),
              lower = matrix(NA,nrow=nscen,ncol=nci),
              upper = matrix(NA,nrow=nscen,ncol=nci))
  for (i in 1:nscen) {
    simmu <- b%*%x[i,]
    # simy <- 1/(1+ exp(-simmu)) # use cloglog instead as below:
    simy <- 1-exp(-exp(simmu))
    res$pe[i] <- mean(simy)
    for (k in 1:nci) {
      cint <- quantile(simy,probs=c( (1-ci[k])/2, (1-(1-ci[k])/2) ) )
      res$lower[i,k] <- cint[1]
      res$upper[i,k] <- cint[2]
    }
  }
  res$lower <- drop(res$lower)
  res$upper <- drop(res$upper)
  res
}

# Simulate first difference of expected probabilities for logit
cloglogsimfd <- function(x,b,ci=0.95,constant=1,xpre=NULL) {
  if (any(class(x)=="counterfactual")&&!is.null(x$model)) {
    xpre <- model.matrix(x$model,x$xpre)     
    x <- model.matrix(x$model,x$x)
  } else {
    if (any(class(x)=="list")) x <- x$x
    if (any(class(x)=="list")) xpre <- x$xpre
    if (is.data.frame(x)) x <- as.matrix(x)
    if (is.data.frame(xpre)) x <- as.matrix(xpre)
    if (!is.matrix(x)) {
      if (is.matrix(b)) {
        x <- t(x)
        if (!is.na(constant)) {
          x <- append(x,1,constant-1)
        }
      } else {
        x <- as.matrix(x)
        if (!is.na(constant)) {
          x <- appendmatrix(x,rep(1,nrow(x)),constant)
        }
      }
    } else {
      if (!is.na(constant)) {
        x <- appendmatrix(x,rep(1,nrow(x)),constant)
      }
    }
    if (!is.matrix(xpre)) {
      if (is.matrix(b)) {
        xpre <- t(xpre)
        if (!is.na(constant)) {
          xpre <- append(xpre,1,constant-1)
        }
      } else {
        xpre <- as.matrix(xpre)
        if (!is.na(constant)) {
          xpre <- appendmatrix(xpre,rep(1,nrow(xpre)),constant)
        }
      }
    } else {
      if (!is.na(constant)) {
        xpre <- appendmatrix(xpre,rep(1,nrow(xpre)),constant)
      }
    }
  }
  
  esims <- nrow(as.matrix(b))
  
  nscen <- nrow(x)
  nci <- length(ci)
  res <- list(pe = rep(NA, nscen),
              lower = matrix(NA,nrow=nscen,ncol=nci),
              upper = matrix(NA,nrow=nscen,ncol=nci))
  for (i in 1:nscen) {
    simmu1 <- b%*%xpre[i,]
    simmu2 <- b%*%x[i,]
    # simy <- 1/(1+ exp(-simmu2)) - 1/(1+ exp(-simmu1)) # use cloglog instead as below:
    simy <- 1-exp(-exp(simmu2)) - (1-exp(-exp(simmu1)))
    res$pe[i] <- mean(simy)
    for (k in 1:nci) {
      cint <- quantile(simy,probs=c( (1-ci[k])/2, (1-(1-ci[k])/2) ) )
      res$lower[i,k] <- cint[1]
      res$upper[i,k] <- cint[2]
    }
  }
  res$lower <- drop(res$lower)
  res$upper <- drop(res$upper)
  res
}



# Simulate relative risk of expected probabilities for logit
cloglogsimrr <- function(x,b,ci=0.95,constant=1,xpre=NULL) {
  if (any(class(x)=="counterfactual")&&!is.null(x$model)) {
    xpre <- model.matrix(x$model,x$xpre)   
    x <- model.matrix(x$model,x$x)
  } else {
    if (any(class(x)=="list")) x <- x$x
    if (any(class(x)=="list")) xpre <- x$xpre
    if (is.data.frame(x)) x <- as.matrix(x)
    if (is.data.frame(xpre)) x <- as.matrix(xpre)
    if (!is.matrix(x)) {
      if (is.matrix(b)) {
        x <- t(x)
        if (!is.na(constant)) {
          x <- append(x,1,constant-1)
        }
      } else {
        x <- as.matrix(x)
        if (!is.na(constant)) {
          x <- appendmatrix(x,rep(1,nrow(x)),constant)
        }
      }
    } else {
      if (!is.na(constant)) {
        x <- appendmatrix(x,rep(1,nrow(x)),constant)
      }
    }
    if (!is.matrix(xpre)) {
      if (is.matrix(b)) {
        xpre <- t(xpre)
        if (!is.na(constant)) {
          xpre <- append(xpre,1,constant-1)
        }
      } else {
        xpre <- as.matrix(xpre)
        if (!is.na(constant)) {
          xpre <- appendmatrix(xpre,rep(1,nrow(xpre)),constant)
        }
      }
    } else {
      if (!is.na(constant)) {
        xpre <- appendmatrix(xpre,rep(1,nrow(xpre)),constant)
      }
    }
  }
  
  esims <- nrow(as.matrix(b))
  
  nscen <- nrow(x)
  nci <- length(ci)
  res <- list(pe = rep(NA, nscen),
              lower = matrix(NA,nrow=nscen,ncol=nci),
              upper = matrix(NA,nrow=nscen,ncol=nci))
  for (i in 1:nscen) {
    simmu1 <- b%*%xpre[i,]
    simmu2 <- b%*%x[i,]
    # simy <- (1/(1+ exp(-simmu2))) / (1/(1+ exp(-simmu1))) # use cloglog instead as below:
    simy <- (1-exp(-exp(simmu2))) / (1-exp(-exp(simmu1)))
    res$pe[i] <- mean(simy)
    for (k in 1:nci) {
      cint <- quantile(simy,probs=c( (1-ci[k])/2, (1-(1-ci[k])/2) ) )
      res$lower[i,k] <- cint[1]
      res$upper[i,k] <- cint[2]
    }
  }
  res$lower <- drop(res$lower)
  res$upper <- drop(res$upper)
  res
}
